#' Plot Gower's distance matrix
#'
#' \code{plot_dist} generates a distance matrix heatmap using the Gower's distance.
#'
#' @param x SummarizedExperiment or DEGdata object,
#' for which differentially expressed/enriched proteins or genes are annotated
#' (output from \code{\link{test_diff}()} and \code{\link{add_rejections}()}).
#' @param significant Logical(1),
#' Whether or not to filter for significant proteins.
#' @param pal Character(1),
#' Sets the color panel (from \pkg{RColorBrewer}).
#' @param pal_rev Logical(1),
#' Whether or not to invert the color palette.
#' @param indicate Character,
#' Sets additional annotation on the top of the heatmap
#' based on columns from the experimental design (colData).
#' @param font_size Integer(1),
#' Sets the size of the labels.
#' @param plot Logical(1),
#' If \code{TRUE} (default) the distance matrix plot is produced.
#' Otherwise (if \code{FALSE}), the data which the
#' distance matrix plot is based on are returned.
#' @param add_values Logical(1).
#' Display the distance value on heatmap.
#' @param value_size Numeric(1).
#' Font size of distance.
#' @param digits Integer(1),
#' indicating the number of decimal places of distance.
#' @param ... Additional arguments for Heatmap function as depicted in
#' \code{\link[ComplexHeatmap]{Heatmap}}
#' @return A heatmap plot (generated by \code{\link[ComplexHeatmap]{Heatmap}})
#'
#' @examples
#' # Load example
#' data(Silicosis_pg)
#' data <- Silicosis_pg
#' data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")
#'
#' # Make SummarizedExperiment
#' ecols <- grep("LFQ.", colnames(data_unique))
#'
#'
#' ## Load experiement design
#' data(Silicosis_ExpDesign)
#' exp_design <- Silicosis_ExpDesign
#' se <- make_se(data_unique, ecols, exp_design)
#'
#' # Filter and normalize
#' filt <- filter_se(se, thr = 0, fraction = 0.4, filter_formula = ~ Reverse != "+" & Potential.contaminant!="+")
#'
#' norm <- normalize_vsn(filt)
#'
#' # Impute missing values using different functions
#' imputed <- impute(norm, fun = "MinProb", q = 0.05)
#'
#' # Test for differentially expressed proteins
#' diff <- test_diff(imputed, "control", "PBS")
#' dep <- add_rejections(diff, alpha = 0.05, lfc = 1)
#'
#' # Plot correlation matrix
#' plot_dist(dep)
#' @export
plot_dist <- function (x, significant = TRUE, pal = "YlOrRd", pal_rev = TRUE,
                       indicate = NULL, font_size = 12, plot = TRUE, add_values = FALSE, value_size = 10, digits = 2, ...)
{
  assertthat::assert_that(inherits(x, "SummarizedExperiment") | inherits(x,"DESeqDataSet"),
                          is.logical(significant), length(significant) == 1, is.character(pal),
                          length(pal) == 1, is.logical(pal_rev), length(pal_rev) ==
                            1, is.numeric(font_size), length(font_size) == 1,
                          is.logical(plot), length(plot) == 1)
  pals <- RColorBrewer::brewer.pal.info %>% rownames_to_column() %>%
    filter(category != "qual")
  if (!pal %in% pals$rowname) {
    stop("'", pal, "' is not a valid color panel",
         " (qualitative panels also not allowed)\n",
         "Run plot_dist() with one of the following 'pal' options: ",
         paste(pals$rowname, collapse = "', '"), "'",
         call. = FALSE)
  }
  if (significant) {
    if(class(x) == "SummarizedExperiment"){
      if (!"significant" %in% colnames(rowData(x, use.names = FALSE))) {
        stop("'significant' column is not present in '",
             deparse(substitute(x)), "'\nRun add_rejections() to obtain the required column",
             call. = FALSE)
      }
      x <- x[rowData(x, use.names = FALSE)$significant,]
      data = assay(x)
    }else if(class(x) == "DEGdata"){
      if (!"significant" %in% colnames(x@test_result)) {
        stop("'significant' column is not present in '",
             deparse(substitute(x)), "'\nRun add_rejections() to obtain the required column",
             call. = FALSE)
      }
      x <- x[x@test_result$significant,]
      data = x@rlg
    }else{
      stop("x must is a DEP or DEG result when significant is TRUE")
    }
  }

  coldata = colData(x)
  plot_dist_util(data, coldata, pal, pal_rev,
                 indicate, font_size, plot, add_values, value_size, digits, ...)


}

## the plot function of distance
plot_dist_util <- function (data, coldata, pal = "YlOrRd", pal_rev = TRUE,
                          indicate = NULL, font_size = 12, plot = TRUE, add_values = FALSE, value_size = 10, digits = 2, ...)
{
  # assertthat::assert_that(
  #   is.character(pal),
  #   length(pal) == 1, is.logical(pal_rev), length(pal_rev) ==
  #     1, is.numeric(font_size), length(font_size) == 1,
  #   is.logical(plot), length(plot) == 1)
  pals <- RColorBrewer::brewer.pal.info %>% rownames_to_column() %>%
    filter(category != "qual")
  # if (!pal %in% pals$rowname) {
  #   stop("'", pal, "' is not a valid color panel",
  #        " (qualitative panels also not allowed)\n",
  #        "Run plot_my_dist() with one of the following 'pal' options: ",
  #        paste(pals$rowname, collapse = "', '"), "'",
  #        call. = FALSE)
  # }
  if (!is.null(indicate)) {
    assertthat::assert_that(is.character(indicate))
    col_data <- coldata
    columns <- colnames(col_data)
    if (any(!indicate %in% columns)) {
      stop("'", paste0(indicate, collapse = "' and/or '"),
           "' column(s) is/are not present in ", deparse(substitute(dep)),
           ".\nValid columns are: '", paste(columns,
                                            collapse = "', '"), "'.", call. = FALSE)
    }
    anno <- coldata %>% dplyr::select(indicate)
    names <- colnames(anno)
    anno_col <- vector(mode = "list", length = length(names))
    names(anno_col) <- names
    for (i in names) {
      var = anno[[i]] %>% unique() %>% sort()
      if (length(var) == 1)
        cols <- c("black")
      if (length(var) == 2)
        cols <- c("orangered", "cornflowerblue")
      if (length(var) < 7 & length(var) > 2)
        cols <- RColorBrewer::brewer.pal(length(var),
                                         "Pastel1")
      if (length(var) >= 7)
        cols <- RColorBrewer::brewer.pal(length(var),
                                         "Set3")
      names(cols) <- var
      anno_col[[i]] <- cols
    }
    ha1 = HeatmapAnnotation(df = anno, col = anno_col, show_annotation_name = TRUE)
  }
  else {
    ha1 <- NULL
  }

  dist_mat <- cluster::daisy(t(data), metric = "gower") %>%
    as.matrix()
  max <- max(dist_mat)
  ht1 = Heatmap(dist_mat, col = circlize::colorRamp2(seq(0, max, ((max)/7)), if (pal_rev) {rev(RColorBrewer::brewer.pal(8, pal))}
                                                     else {RColorBrewer::brewer.pal(8, pal)
                                                     }), heatmap_legend_param = list(color_bar = "continuous", legend_direction = "horizontal", legend_width = unit(5, "cm"), title_position = "topcenter"), name = "Gower's distance", column_names_gp = gpar(fontsize = font_size),
                row_names_gp = gpar(fontsize = font_size), top_annotation = ha1, cell_fun = if(add_values) {function(j, i, x, y, width, height, fill) {
                  grid.text(sprintf(paste("%.", digits, "f", sep = ""), dist_mat[i, j]), x, y, gp = gpar(fontsize = value_size))}} else {NULL},
                ...)
  if (plot) {
    draw(ht1, heatmap_legend_side = "top")
  }
  else {
    df <- as.data.frame(dist_mat)
    return(df)
  }
}


#' Plot correlation matrix
#'
#' \code{plot_cor} generates a Pearson correlation matrix.
#'
#' @param x SummarizedExperiment,
#' Data object for which differentially enriched proteins are annotated
#' (output from \code{\link{test_diff}()} and \code{\link{add_rejections}()}).
#' @param significant Logical(1),
#' Whether or not to filter for significant proteins.
#' @param lower Integer(1),
#' Sets the lower limit of the color scale.
#' @param upper Integer(1),
#' Sets the upper limit of the color scale.
#' @param pal Character(1),
#' Sets the color panel (from \pkg{RColorBrewer}).
#' @param pal_rev Logical(1),
#' Whether or not to invert the color palette.
#' @param indicate Character,
#' Sets additional annotation on the top of the heatmap
#' based on columns from the experimental design (colData).
#' @param font_size Integer(1),
#' Sets the size of the labels.
#' @param plot Logical(1),
#' If \code{TRUE} (default) the correlation matrix plot is produced.
#' Otherwise (if \code{FALSE}), the data which the
#' correlation matrix plot is based on are returned.
#' @param add_values Logical(1).
#' Display the distance value on heatmap.
#' @param value_size Numeric(1).
#' Font size of distance.
#' @param digits Integer(1),
#' indicating the number of decimal places of distance.
#' @param ... Additional arguments for Heatmap function as depicted in
#' \code{\link[ComplexHeatmap]{Heatmap}}
#' @return A heatmap plot (generated by \code{\link[ComplexHeatmap]{Heatmap}})
#'
#' @examples
#' # Load example
#' data(Silicosis_pg)
#' data <- Silicosis_pg
#' data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")
#'
#' # Make SummarizedExperiment
#' ecols <- grep("LFQ.", colnames(data_unique))
#'
#'
#' ## Load experiement design
#' data(Silicosis_ExpDesign)
#' exp_design <- Silicosis_ExpDesign
#' se <- make_se(data_unique, ecols, exp_design)
#'
#' # Filter and normalize
#' filt <- filter_se(se, thr = 0, fraction = 0.4, filter_formula = ~ Reverse != "+" & Potential.contaminant!="+")
#'
#' norm <- normalize_vsn(filt)
#'
#' # Impute missing values using different functions
#' imputed <- impute(norm, fun = "MinProb", q = 0.05)
#'
#' # Test for differentially expressed proteins
#' diff <- test_diff(imputed, "control", "PBS")
#' dep <- add_rejections(diff, alpha = 0.05, lfc = 1)
#'
#' # Plot correlation matrix
#' plot_cor(dep)
#' @export
plot_cor <- function (x, significant = TRUE, lower = -1, upper = 1, pal = "PRGn",
                      pal_rev = FALSE, indicate = NULL, font_size = 12, plot = TRUE,  add_values = FALSE, value_size = 10, digits = 2,
                      ...)
{
  assertthat::assert_that(inherits(x, "SummarizedExperiment"),
                          is.logical(significant), length(significant) == 1, is.numeric(lower),
                          length(lower) == 1, is.numeric(upper), length(upper) ==
                            1, is.character(pal), length(pal) == 1, is.logical(pal_rev),
                          length(pal_rev) == 1, is.numeric(font_size), length(font_size) ==
                            1, is.logical(plot), length(plot) == 1)
  if (!(lower >= -1 & upper >= -1 & lower <= 1 & upper <= 1)) {
    stop("'lower' and/or 'upper' arguments are not valid\n         Run plot_cor() with 'lower' and 'upper' between -1 and 1",
         call. = FALSE)
  }
  pals <- RColorBrewer::brewer.pal.info %>% rownames_to_column() %>%
    filter(category != "qual")
  if (!pal %in% pals$rowname) {
    stop("'", pal, "' is not a valid color panel",
         " (qualitative panels also not allowed)\n",
         "Run plot_cor() with one of the following 'pal' options: ",
         paste(pals$rowname, collapse = "', '"), "'",
         call. = FALSE)
  }
  if (any(is.na(assay(x)))) {
    stop("Missing values in '", deparse(substitute(x)),
         "'. Use plot_dist() instead")
  }

  if (significant) {
    if(class(x) == "SummarizedExperiment"){
      if (!"significant" %in% colnames(rowData(x, use.names = FALSE))) {
        stop("'significant' column is not present in '",
             deparse(substitute(x)), "'\nRun add_rejections() to obtain the required column",
             call. = FALSE)
      }
      x <- x[rowData(x, use.names = FALSE)$significant,
      ]
    }else if(class(x) == "DEGdata"){
      if (!"significant" %in% colnames(x@test_result)) {
        stop("'significant' column is not present in '",
             deparse(substitute(x)), "'\nRun add_rejections() to obtain the required column",
             call. = FALSE)
      }
      x <- x[x@test_result$significant,]
    }else{
      stop("x must is a DEP or DEG result when significant is TRUE")
    }
  }
  data = assay(x)
  coldata = colData(x)

  plot_cor_util(data, coldata, lower, upper, pal,
                pal_rev, indicate, font_size, plot, add_values, value_size, digits,
                ...)
}

plot_cor_util <- function (data, coldata, lower = -1, upper = 1, pal = "PRGn",
                         pal_rev = FALSE, indicate = NULL, font_size = 12, plot = TRUE, add_values = FALSE, value_size = 10, digits = 2,
                         ...)
{
  # assertthat::assert_that(
  #   is.numeric(lower),
  #   length(lower) == 1, is.numeric(upper), length(upper) ==
  #     1, is.character(pal), length(pal) == 1, is.logical(pal_rev),
  #   length(pal_rev) == 1, is.numeric(font_size), length(font_size) ==
  #     1, is.logical(plot), length(plot) == 1)
  # if (!(lower >= -1 & upper >= -1 & lower <= 1 & upper <= 1)) {
  #   stop("'lower' and/or 'upper' arguments are not valid\n         Run plot_cor() with 'lower' and 'upper' between -1 and 1",
  #        call. = FALSE)
  # }
  # pals <- RColorBrewer::brewer.pal.info %>% rownames_to_column() %>%
  #   filter(category != "qual")
  # if (!pal %in% pals$rowname) {
  #   stop("'", pal, "' is not a valid color panel",
  #        " (qualitative panels also not allowed)\n",
  #        "Run plot_cor() with one of the following 'pal' options: ",
  #        paste(pals$rowname, collapse = "', '"), "'",
  #        call. = FALSE)
  # }
  # if (any(is.na(data))) {
  #   stop("Missing values in '", deparse(substitute(data)),
  #        "'. Use plot_dist() instead")
  # }
  if (!is.null(indicate)) {
    assertthat::assert_that(is.character(indicate))
    col_data <- coldata %>% as.data.frame()
    columns <- colnames(col_data)
    if (any(!indicate %in% columns)) {
      stop("'", paste0(indicate, collapse = "' and/or '"),
           "' column(s) is/are not present in ", deparse(substitute(dep)),
           ".\nValid columns are: '", paste(columns,
                                            collapse = "', '"), "'.", call. = FALSE)
    }
    anno <- coldata %>% data.frame() %>% dplyr::select(indicate)
    names <- colnames(anno)
    anno_col <- vector(mode = "list", length = length(names))
    names(anno_col) <- names
    for (i in names) {
      var = anno[[i]] %>% unique() %>% sort()
      if (length(var) == 1)
        cols <- c("black")
      if (length(var) == 2)
        cols <- c("orangered", "cornflowerblue")
      if (length(var) < 7 & length(var) > 2)
        cols <- RColorBrewer::brewer.pal(length(var),
                                         "Pastel1")
      if (length(var) >= 7)
        cols <- RColorBrewer::brewer.pal(length(var),
                                         "Set3")
      names(cols) <- var
      anno_col[[i]] <- cols
    }
    ha1 = HeatmapAnnotation(df = anno, col = anno_col, show_annotation_name = TRUE)
  }
  else {
    ha1 <- NULL
  }
  # if (significant) {
  #     if (!"significant" %in% colnames(rowData(dep, use.names = FALSE))) {
  #         stop("'significant' column is not present in '",
  #             deparse(substitute(dep)), "'\nRun add_rejections() to obtain the required column",
  #             call. = FALSE)
  #     }
  #     data <- dep[data$significant,
  #         ]
  # }
  cor_mat <- cor(data)
  ht1 = Heatmap(cor_mat, col = circlize::colorRamp2(seq(lower,
                                                        upper, ((upper - lower)/7)), if (pal_rev) {
                                                          rev(RColorBrewer::brewer.pal(8, pal))
                                                        }
                                                    else {
                                                      RColorBrewer::brewer.pal(8, pal)
                                                    }), heatmap_legend_param = list(color_bar = "continuous",
                                                                                    legend_direction = "horizontal", legend_width = unit(5,
                                                                                                                                         "cm"), title_position = "topcenter"),
                name = "Pearson correlation", column_names_gp = gpar(fontsize = font_size),
                row_names_gp = gpar(fontsize = font_size), top_annotation = ha1, cell_fun = if(add_values) {function(j, i, x, y, width, height, fill) {
                  grid.text(sprintf(paste("%.", digits, "f", sep = ""), cor_mat[i, j]), x, y, gp = gpar(fontsize = value_size))}} else {NULL},
                ...)
  if (plot) {
    draw(ht1, heatmap_legend_side = "top")
  }
  else {
    df <- as.data.frame(cor_mat)
    return(df)
  }
}





#' Plot PCA
#'
#' \code{plot_pca} generates a PCA plot using the top variable features
#'
#' @param object SummarizedExperiment (or DEGdata) object,
#' Data object for which differentially enriched proteins are annotated
#' (output from \code{\link{test_diff}()} (or \code{\link{test_diff_deg}()}) and \code{\link{add_rejections}()}).
#' @param x Integer(1),
#' Sets the principle component to plot on the x-axis.
#' @param y Integer(1),
#' Sets the principle component to plot on the y-axis.
#' @param indicate Character,
#' Sets the color, shape and facet_wrap of the plot
#' based on columns from the experimental design (colData).
#' @param label Logical,
#' Whether or not to add sample labels.
#' @param n Integer(1),
#' Sets the number of top variable proteins to consider.
#' @param point_size Integer(1),
#' Sets the size of the points.
#' @param label_size Integer(1),
#' Sets the size of the labels.
#' @param plot Logical(1),
#' If \code{TRUE} (default) the PCA plot is produced.
#' Otherwise (if \code{FALSE}), the data which the
#' PCA plot is based on are returned.
#' @param features Character(1), the feature name in plot title, could be "proteins","genes", default is "features"
#' @param if_square Logical(1), if TRUE plot in a
#' @return A scatter plot (generated by \code{\link[ggplot2]{ggplot}}).
#'
#' @examples
#' # Load example
#' data(Silicosis_pg)
#' data <- Silicosis_pg
#' data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")
#'
#' # Differential test
#' ecols <- grep("LFQ.", colnames(data_unique))
#' se <- make_se_parse(data_unique, ecols,mode = "delim")
#' filt <- filter_se(se, thr = 0, fraction = 0.4, filter_formula = ~ Reverse != "+" & Potential.contaminant!="+")
#' norm <- normalize_vsn(filt)
#' imputed <- impute(norm, fun = "MinProb", q = 0.05)
#'
#' # UMAP plot
#' plot_pca(imputed)
#' @export
plot_pca <- function(object, x = 1, y = 2, indicate = c("condition", "replicate"), label = FALSE, n = 500, point_size = 4,
                     label_size = 3, plot = TRUE, features = "features", if_square = FALSE){

  if (is.integer(x))
    x <- as.numeric(x)
  if (is.integer(y))
    y <- as.numeric(y)
  if (is.integer(n))
    n <- as.numeric(n)
  if (is.integer(point_size))
    point_size <- as.numeric(point_size)
  if (is.integer(label_size))
    label_size <- as.numeric(label_size)

  assertthat::assert_that(
    is.numeric(x), length(x) == 1, is.numeric(y), length(y) ==
      1, is.numeric(n), length(n) == 1, is.character(indicate),
    is.logical(label), is.numeric(point_size), length(point_size) ==
      1, is.numeric(label_size), length(label_size) ==
      1, is.logical(plot), length(plot) == 1,
    is.character(features), length(features) == 1,
    is.logical(if_square), length(if_square) == 1)

  data = assay(object)
  coldata = colData(object)
  if (x > ncol(data) | y > ncol(data)) {
    stop(paste0("'x' and/or 'y' arguments are not valid\n",
                "Run plot_pca() with 'x' and 'y' <= ", ncol(data),
                "."), call. = FALSE)
  }
  if (n > nrow(data)) {
    stop(paste0("'n' argument is not valid.\n", "Run plot_pca() with 'n' <= ",
                nrow(data), "."), call. = FALSE)
  }

  # for shiny indicate can select the third indicate
  mycolData <- coldata
  mycolData$Condition = mycolData$condition
  mycolData$Replicate = mycolData$replicate

  columns <- colnames(mycolData)
  if (!is.null(indicate)) {
    if (length(indicate) > 3) {
      stop("Too many features in 'indicate'\n        Run plot_pca() with a maximum of 3 indicate features")
    }
    if (any(!indicate %in% columns)) {
      stop(paste0("'", paste0(indicate, collapse = "' and/or '"),
                  "' column(s) is/are not present in ", deparse(substitute(data)),
                  ".\nValid columns are: '", paste(columns,
                                                   collapse = "', '"), "'."), call. = FALSE)
    }
  }
  var <- apply(data, 1, sd)
  df <- data[order(var, decreasing = TRUE)[seq_len(n)],]

  pca <- prcomp(t(df), scale = FALSE)
  pca_df <- pca$x %>% data.frame() %>% rownames_to_column() %>%
    dplyr::left_join(., data.frame(mycolData), by = c(rowname = "ID"))
  percent <- round(100 * pca$sdev^2/sum(pca$sdev^2), 1)
  for (feat in indicate) {
    pca_df[[feat]] <- as.factor(pca_df[[feat]])
  }

  limit = unlist(pca_df[,c(paste0("PC", x), paste0("PC", y))])
  p <- ggplot(pca_df, aes(get(paste0("PC", x)), get(paste0("PC", y)))) +
    labs(title = paste0("PCA plot - top ", n, " variable ", features), x = paste0("PC", x, ": ", percent[x], "%"), y = paste0("PC", y, ": ", percent[y], "%")) +
    coord_fixed() +
    theme_DEP1()
  if(if_square){
    p <- p +
      scale_x_continuous(limits = range(limit)) +
      scale_y_continuous(limits = range(limit))
  }

  if (length(indicate) == 0) {
    p <- p + geom_point(size = point_size)
  }
  if (length(indicate) == 1) {
    p <- p + geom_point(aes(col = pca_df[[indicate[1]]]),
                        size = point_size) + labs(col = indicate[1])
  }
  if (length(indicate) == 2) {
    p <- p + geom_point(aes(col = pca_df[[indicate[1]]],
                            shape = pca_df[[indicate[2]]]), size = point_size) +
      labs(col = indicate[1], shape = indicate[2])

    if(length(unique(pca_df[[indicate[2]]])) > 6){
      p <- p + scale_shape_manual( values = seq(1, length(unique(pca_df[[indicate[2]]]))) )
    }
  }
  if (length(indicate) == 3) {
    p <- p + geom_point(aes(col = pca_df[[indicate[1]]],
                            shape = pca_df[[indicate[2]]]), size = point_size) +
      facet_wrap(~pca_df[[indicate[3]]]) +
      labs(col = indicate[1], shape = indicate[2])

    if(length(unique(pca_df[[indicate[2]]])) > 6){
      p <- p + scale_shape_manual( values = seq(1, length(unique(pca_df[[indicate[2]]]))) )
    }
  }
  if (label) {
    p <- p + geom_text(aes(label = rowname), size = label_size)
  }
  if (plot) {
    return(p)
  }
  else {
    df <- pca_df %>% dplyr::select(rowname, paste0("PC", c(x,
                                                           y)), match(indicate, colnames(pca_df)))
    colnames(df)[1] <- "sample"
    return(df)
  }

}


#' Plot t-Sne
#'
#' \code{plot_Tsne} generates a t-Sne plot using the global variable features
#'
#' @param perplexity 	Perplexity parameter (should not be bigger than 3 * perplexity < nrow(X) - 1. See \link[Rtsne]{Rtsne}().
#' @param theseed integer(1),Set random seed.
#' @param ... Other parameters to \link[Rtsne]{Rtsne}().
#' @inheritParams plot_pca
#'
#' @return A scatter plot (generated by \code{\link[ggplot2]{ggplot}}).
#'
#' @importFrom Rtsne Rtsne
#' @export
#'
#' @examples
#' # Load example
#' data(Silicosis_pg)
#' data <- Silicosis_pg
#' data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")
#'
#' # Differential test
#' ecols <- grep("LFQ.", colnames(data_unique))
#' se <- make_se_parse(data_unique, ecols,mode = "delim")
#' filt <- filter_se(se, thr = 0, fraction = 0.4, filter_formula = ~ Reverse != "+" & Potential.contaminant!="+")
#' norm <- normalize_vsn(filt)
#' imputed <- impute(norm, fun = "MinProb", q = 0.05)
#'
#' # t-SNE plot
#' plot_Tsne(imputed)
#'
plot_Tsne <- function (object, indicate = c("condition",
                                         "replicate"), label = FALSE,
                       n = nrow(object),
                       point_size = 4, label_size = 3, plot = TRUE, if_square = FALSE, features = "features",
                       perplexity = 2, theseed = NULL,...)
{
  if (is.numeric(n))
    n <- as.integer(n)
  if (is.numeric(theseed))
    theseed <- as.integer(theseed)
  if (is.integer(point_size))
    point_size <- as.numeric(point_size)
  if (is.integer(label_size))
    label_size <- as.numeric(label_size)
  assertthat::assert_that(is.numeric(n), length(n) == 1,
                          is.character(indicate),
                          is.logical(label), is.numeric(point_size), length(point_size) ==1,
                          is.numeric(label_size), length(label_size) ==1,
                          is.logical(plot), length(plot) == 1,
                          is.numeric(perplexity),length(perplexity)==1,
                          is.null(theseed)|is.integer(theseed))

  coldata <- colData(object)
  data <- assay(object)

  if (n > nrow(data)) {
    stop(paste0("'n' argument is not valid.\n", "Run plot_Tsne_for_rnaseq() with 'n' <= ",
                nrow(data), "."), call. = FALSE)
  }

  mycolData <- coldata
  mycolData$Condition = mycolData$condition
  mycolData$Replicate = mycolData$replicate
  columns <- colnames(mycolData)

  if (!is.null(indicate)) {
    if (length(indicate) > 3) {
      stop("Too many features in 'indicate'\n        Run plot_Tsne_for_rnaseq() with a maximum of 3 indicate features")
    }
    if (any(!indicate %in% columns)) {
      stop(paste0("'", paste0(indicate, collapse = "' and/or '"),
                  "' column(s) is/are not present in ", deparse(substitute(data)),
                  ".\nValid columns are: '", paste(columns,
                                                   collapse = "', '"), "'."), call. = FALSE)
    }
  }
  var <- apply(data, 1, sd)
  n = min(n,nrow(data))
  # n_neighbors = min( nrow(t(df)),n_neighbors)
  df <- data[order(var, decreasing = TRUE)[seq_len(n)], ]

  if(is.null(theseed)) {
    theseed = sample(1:2^15,1)
    message("The random seed of umap is ", theseed)
  }
  set.seed(theseed)

  Tsne_res = Rtsne::Rtsne(t(df),perplexity = perplexity)
  Tsne_df  = Tsne_res$Y %>% data.frame() %>% `rownames<-`(colnames(df)) %>%rownames_to_column() %>%
    left_join(., data.frame(mycolData), by = c(rowname = "ID"))
  colnames(Tsne_df)[2:3] = c("x","y")

  limit = unlist(Tsne_df[,2:3])
  p = ggplot(Tsne_df,aes(x,y)) + labs(title = paste0("T-SNE plot top ",n, " variable ",features))+
    coord_fixed() +
    theme_DEP1()
  if(if_square){
    p <- p +
      scale_x_continuous(limits = range(limit)) +
      scale_y_continuous(limits = range(limit))
  }

  if (length(indicate) == 0) {
    p <- p + geom_point(size = point_size)
  }
  if (length(indicate) == 1) {
    p <- p + geom_point(aes(col = as.factor(Tsne_df[[indicate[1]]])),
                        size = point_size) + labs(col = indicate[1])
  }
  if (length(indicate) == 2) {
    p <- p + geom_point(aes(col = as.factor(Tsne_df[[indicate[1]]]),
                            shape = as.factor(Tsne_df[[indicate[2]]])), size = point_size) +
      labs(col = indicate[1], shape = indicate[2])

    if(length(unique(Tsne_df[[indicate[2]]])) > 6){
      p <- p + scale_shape_manual( values = seq(1, length(unique(Tsne_df[[indicate[2]]]))) )
    }
  }
  if (length(indicate) == 3) {
    p <- p + geom_point(aes(col = as.factor(Tsne_df[[indicate[1]]]),
                            shape = as.factor(Tsne_df[[indicate[2]]])), size = point_size) +
      facet_wrap(~Tsne_df[[indicate[3]]]) +
      labs(col = indicate[1], shape = indicate[2])

    if(length(unique(Tsne_df[[indicate[2]]])) > 6){
      p <- p + scale_shape_manual( values = seq(1, length(unique(Tsne_df[[indicate[2]]]))) )
    }
  }
  if (label) {
    p <- p + geom_text(aes(label = rowname), size = label_size)
  }
  if (plot) {
    p <- p + labs(x = "", y = "")
    return(p)
  }
  else {
    df <- Tsne_df %>% dplyr::select(rowname, c("x","y"), match(indicate, colnames(Tsne_df)))
    colnames(df)[1] <- "sample"
    return(df)
  }

}


#indicate can be one , two, three element, the third means facet_wrap
# if_square: if xlim == ylim, default FALSE
#' Plot Tsne
#'
#' \code{plot_umap} generates a umap plot using global variable features though \link[umap]{umap}.
#'
#' @param n_neighbors Integer(1), constraining the size of the local neighborhood UMAP will look at.
#' @param theseed Integer(1),Set random_state i \link[umap]{umap}.
#' #' @param ...
#' @inheritParams plot_pca
#'
#' @return A scatter plot (generated by \code{\link[ggplot2]{ggplot}}).
#'
#' @importFrom umap umap
#' @export
#'
#' @examples
#' # Load example
#' data(Silicosis_pg)
#' data <- Silicosis_pg
#' data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")
#'
#' # Differential test
#' ecols <- grep("LFQ.", colnames(data_unique))
#' se <- make_se_parse(data_unique, ecols,mode = "delim")
#' filt <- filter_se(se, thr = 0, fraction = 0.4, filter_formula = ~ Reverse != "+" & Potential.contaminant!="+")
#' norm <- normalize_vsn(filt)
#' imputed <- impute(norm, fun = "MinProb", q = 0.05)
#'
#' # UMAP plot
#' plot_umap(imputed)
plot_umap <- function (object, indicate = c("condition",
                                         "replicate"), label = FALSE, n = nrow(object),
                       point_size = 4, label_size = 3, plot = TRUE, if_square = FALSE, features = "features",
                       n_neighbors=3, theseed=NULL, ...){
  if (is.numeric(n))
    n <- as.integer(n)
  if (is.numeric(n_neighbors))
    n_neighbors <- as.integer(n_neighbors)
  if (is.numeric(theseed))
    theseed <- as.integer(theseed)
  if (is.integer(point_size))
    point_size <- as.numeric(point_size)
  if (is.integer(label_size))
    label_size <- as.numeric(label_size)
  assertthat::assert_that(is.character(indicate),
                          is.logical(label), is.numeric(point_size), length(point_size) ==1,
                          is.numeric(label_size), length(label_size) ==1,
                          is.logical(plot), length(plot) == 1,
                          is.integer(n_neighbors),length(n_neighbors)==1,
                          is.integer(n),length(n) ==1,
                          is.null(theseed)|is.integer(theseed))

  coldata <- colData(object)
  data <- assay(object)

  if (n > nrow(data)) {
    stop(paste0("'n' argument is not valid.\n", "Run plot_umap_for_rnaseq() with 'n' <= ",
                nrow(data), "."), call. = FALSE)
  }
  if (n_neighbors > ncol(data)){
    stop(paste0("'n_neighbors' argument is not valid.\n", "Run plot_umap_for_rnaseq() with 'n_neighbors' <= ",
                ncol(data), "."), call. = FALSE)
  }
  # for shiny indicate can select the third indicate
  mycolData <- coldata
  mycolData$Condition = mycolData$condition
  mycolData$Replicate = mycolData$replicate
  columns <- colnames(mycolData)

  # cat("vvv")
  if (!is.null(indicate)) {
    if (length(indicate) > 3) {
      stop("Too many features in 'indicate'\n        Run plot_umap_for_rnaseq() with a maximum of 3 indicate features")
    }
    if (any(!indicate %in% columns)) {
      stop(paste0("'", paste0(indicate, collapse = "' and/or '"),
                  "' column(s) is/are not present in ", deparse(substitute(data)),
                  ".\nValid columns are: '", paste(columns,
                                                   collapse = "', '"), "'."), call. = FALSE)
    }
  }
  # cat("ddd")
  var <- apply(data, 1, sd)
  n = min(n,nrow(data))
  df <- data[order(var, decreasing = TRUE)[seq_len(n)], ]

  if(is.null(theseed)) {
    theseed = sample(1:2^15,1)
    message("The random_state of umap is ", theseed)
  }

  n_neighbors = min( nrow(t(df)),n_neighbors)
  umap_res = umap::umap(t(df),n_neighbors= n_neighbors, random_state=theseed, ...)
  umap_df  = umap_res$layout %>% data.frame() %>% rownames_to_column() %>%
    left_join(., data.frame(mycolData), by = c(rowname = "ID"))
  colnames(umap_df)[2:3] = c("x","y")

  limit = unlist(umap_df[,2:3])
  # cat("eee")
  p = ggplot(umap_df,aes(x,y)) + labs(title = paste0("UMAP plot top ",n, " variable", features))+
    coord_fixed() +
    theme_DEP1()
  if(if_square){
    p <- p +
      scale_x_continuous(limits = range(limit)) +
      scale_y_continuous(limits = range(limit))
  }

  if (length(indicate) == 0) {
    p <- p + geom_point(size = point_size)
  }
  if (length(indicate) == 1) {
    p <- p + geom_point(aes(col = umap_df[[indicate[1]]]),
                        size = point_size) + labs(col = as.factor(indicate[1]))
  }
  if (length(indicate) == 2) {
    p <- p + geom_point(aes(col = as.factor(umap_df[[indicate[1]]]),
                            shape = as.factor(umap_df[[indicate[2]]])), size = point_size) +
      labs(col = indicate[1], shape = indicate[2])

    if(length(unique(umap_df[[indicate[2]]])) > 6){
      p <- p + scale_shape_manual( values = seq(1, length(unique(umap_df[[indicate[2]]]))) )
    }
  }
  if (length(indicate) == 3) {
    p <- p + geom_point(aes(col = as.factor(umap_df[[indicate[1]]]),
                            shape = as.factor(umap_df[[indicate[2]]])), size = point_size) +
      facet_wrap(~umap_df[[indicate[3]]]) +
      labs(col = indicate[1], shape = indicate[2])

    if(length(unique(umap_df[[indicate[2]]])) > 6){
      p <- p + scale_shape_manual( values = seq(1, length(unique(umap_df[[indicate[2]]]))) )
    }
  }
  if (label) {
    p <- p + geom_text(aes(label = rowname), size = label_size)
  }
  if (plot) {
    p <- p + labs(x = "", y = "")
    return(p)
  }
  else {
    df <- umap_df %>% dplyr::select(rowname, c("x","y"), match(indicate, colnames(umap_df)))
    colnames(df)[1] <- "sample"
    return(df)
  }
}


#' MA-plot of RNA expression data
#'
#' Plot MA-plot for a differential test result out from \code{\link{test_diff}()}.
#'
#' @param x DEGdata object,
#' Test result for which differentially expressed/enriched proteins or genes are annotated
#' (output from \code{\link{test_diff}()}).
#' @param contrast Character.
#' Only plot on certain contrasts.
#' @param point_alpha Numeric(1), between 0 to 1,
#' transparency value for the points (0 = transparent, 1 = opaque)
#' @param sig_color Charecter(1),
#' Color to use to mark differentially expressed genes. Defaults to red
#' @param hlines Numeric(1),
#' The y coordinate (in absolute value) where to draw horizontal lines,
#' optional
#' @param title A title for the plot, optional
#' @param xlab X axis label, defaults to "mean of normalized counts - log10 scale"
#' @param ylim Vector of two numeric values, Y axis limits to restrict the view.
#' @param add_rug Logical, whether to add rug plots in the margins
#' @param intgenes Vector of genes of interest, highlighted on plot by \code{intgenes_color}.
#' Gene symbols if a \code{symbol} column is provided in \code{x},
#' or else the identifiers specified in the row names
#' @param intgenes_color Character(1), The color to use to mark the genes on the main plot.
#' @param intgenes_size The point size of \code{intgenes}.
#' @param labels_intgenes Logical, whether to add the gene identifiers/names close
#' to the marked plots.
#' @param label_fontsize Numeric(1),
#' The fontsize of \code{labels_intgenes}
#' @param labels_repel Logical, whether to use \code{geom_text_repel} for placing the
#' labels on the features to mark
#' @return
#' A ggplot object
#'
#' @export
#\' @importFrom ggrepel geom_label_repel
plot_ma_RNA <- function(x,
                        contrast = get_contrast(x)[1],
                        point_alpha = 0.2,
                        sig_color = "red3",
                        hlines = NULL,
                        title = contrast,
                        xlab = "mean of normalized counts - log10 scale",
                        ylim = NULL,
                        add_rug = FALSE,
                        intgenes = NULL,
                        intgenes_color = "steelblue4",
                        intgenes_size = 3,
                        # labels_intgenes = TRUE,
                        label_fontsize = 3,
                        labels_repel = FALSE) {
  assertthat::assert_that(class(x) == "DEGdata",
                          is.character(contrast), length(contrast) == 1,
                          is.character(sig_color), length(sig_color) == 1,
                          is.character(intgenes_color), length(intgenes_color) == 1,
                          is.character(xlab), is.character(xlab) == 1, is.null(ylim) || is.numeric(ylim),
                          # is.logical(labels_intgenes), length(labels_intgenes) == 1,
                          is.logical(labels_repel), length(labels_repel) == 1,
                          is.numeric(point_alpha), length(point_alpha) == 1,
                          is.numeric(intgenes_size), length(intgenes_size) == 1,
                          is.numeric(label_fontsize), length(label_fontsize) == 1
                          )
  x = x@test_result

  if (length(grep("_significant", colnames(x))) < 1) {
    stop(paste0("'[contrast]_significant' columns are not present in '",
                deparse(substitute(x)), "'.\nRun add_rejections() to obtain the required columns."),
         call. = FALSE)
  }
  if (length(grep(paste("^", contrast, "_diff", sep = ""),
                  colnames(x))) == 0) {
    valid_cntrsts <-  get_contrast(x)
    valid_cntrsts_msg <- paste0("Valid contrasts are: '",
                                paste0(valid_cntrsts, collapse = "', '"), "'")
    stop("Not a valid contrast, please run `plot_ma_RNA()` with a valid contrast as argument\n",
         valid_cntrsts_msg, call. = FALSE)
  }

  ma_df <- data.frame(
    mean = x[ , paste(contrast, "baseMean", sep = "_")],
    lfc = x[ , paste(contrast, "diff", sep = "_")],
    padj = x[ , paste(contrast, "p.adj", sep = "_")],
    # isDE = ifelse(is.na(x[ , paste(contrast, "padj", sep = "_")]), FALSE, x[ , paste(contrast, "padj", sep = "_")] < FDR & abs(x[ , paste(contrast, "log2FoldChange", sep = "_")]) >= log2fc),
    isDE = x[ , paste(contrast, "significant", sep = "_")],
    ID = rownames(x)
  )
  rownames(ma_df) = rownames(x)
  ma_df$name = rownames(ma_df)

  if(!is.null(intgenes) && length(intersect(intgenes, ma_df$name)) > 0 ){
    labels_intgenes = T
    intgenes = intersect(intgenes, ma_df$name)
  }else{
    labels_intgenes = F
  }

  plot_ma_utili(ma_df, point_alpha, sig_color, hlines, title,
                xlab, ylim, add_rug,
                intgenes, intgenes_color, intgenes_size,
                labels_intgenes, labels_repel, label_fontsize)
}


#' MA-plot of quantity data
#'
#' Plot MA-plot for a differential test result.
#'
#' @param x SummarizedExperiment,
#' Data object for which differentially enriched proteins are annotated
#' (output from \code{\link{test_diff}()} and \code{\link{test_diff_deg}()}).
#' @param intgenes Vector of name, highlighted on plot by \code{intgenes_color}
#' column is provided in \code{x}, or else the identifiers specified in the
#' row names
#' @inheritParams plot_ma_RNA
#' @return
#' A ggplot object
#' @export
#'
#\' @importFrom ggrepel geom_label_repel
plot_ma_pro <- function(x,
                        contrast,
                        point_alpha = 0.2,
                        sig_color = "red3",
                        hlines = NULL,
                        title = contrast,
                        xlab = "mean of normalized counts - log10 scale",
                        ylim = NULL,
                        add_rug = FALSE,
                        intgenes = NULL,
                        intgenes_color = "steelblue4",
                        intgenes_size = 3,
                        # labels_intgenes = TRUE,
                        label_fontsize = 3,
                        labels_repel = FALSE) {
  assertthat::assert_that(inherits(x,"SummarizedExperiment"),
                          is.character(contrast), length(contrast) == 1,
                          is.character(sig_color), length(sig_color) == 1,
                          is.character(intgenes_color), length(intgenes_color) == 1,
                          is.character(xlab), length(xlab) == 1,
                          is.null(intgenes) || is.character(intgenes),
                          is.null(ylim) || is.numeric(ylim), is.null(hlines) || is.numeric(hlines),
                          # is.logical(labels_intgenes), length(labels_intgenes) == 1,
                          is.logical(labels_repel), length(labels_repel) == 1,
                          is.numeric(point_alpha), length(point_alpha) == 1,
                          is.numeric(intgenes_size), length(intgenes_size) == 1,
                          is.numeric(label_fontsize), length(label_fontsize) == 1
  )

  dep = x
  x = rowData(dep)
  coldata = colData(dep)
  if (length(grep("_significant", colnames(x))) < 1) {
    stop(paste0("'[contrast]_significant' columns are not present in '",
                deparse(substitute(dep)), "'.\nRun add_rejections() to obtain the required columns."),
         call. = FALSE)
  }
  if (length(grep(paste("^", contrast, "_diff", sep = ""),
                  colnames(x))) == 0) {
    valid_cntrsts <- get_contract(dep)
    valid_cntrsts_msg <- paste0("Valid contrasts are: '",
                                paste0(valid_cntrsts, collapse = "', '"), "'")
    stop("Not a valid contrast, please run `plot_ma_pro()` with a valid contrast as argument\n",
         valid_cntrsts_msg, call. = FALSE)
  }

  theconditions <- strsplit(contrast,"_vs_")[[1]]
  means <- assay(dep)[,which(coldata$condition %in% theconditions)] %>% rowMeans()
  ma_df <- data.frame(
    mean = means,
    lfc = x[ , paste(contrast, "diff", sep = "_")],
    padj = x[ , paste(contrast, "p.adj", sep = "_")],
    # isDE = ifelse(is.na(x[ , paste(contrast, "padj", sep = "_")]), FALSE, x[ , paste(contrast, "padj", sep = "_")] < FDR & abs(x[ , paste(contrast, "log2FoldChange", sep = "_")]) >= log2fc),
    isDE = x[ , paste(contrast, "significant", sep = "_")],
    ID = rownames(x)
  )

  rownames(ma_df) = rownames(x)
  ma_df$name = rownames(ma_df)

  if(!is.null(intgenes) && length(intersect(intgenes, ma_df$name)) > 0 ){
    labels_intgenes = T
    intgenes = intersect(intgenes, ma_df$name)
  }else{
    labels_intgenes = F
  }

  # cat("weg")
  plot_ma_utili(ma_df, point_alpha, sig_color, hlines, title,
                xlab, ylim, add_rug,
                intgenes, intgenes_color, intgenes_size,
                labels_intgenes, labels_repel, label_fontsize)
}


plot_ma_utili <- function(ma_df,
                          point_alpha = 0.2,
                          sig_color = "red3",
                          hlines = NULL,
                          title = contrast,
                          xlab = "mean of normalized counts - log10 scale",
                          ylim = NULL,
                          add_rug = FALSE,
                          intgenes = NULL,
                          intgenes_color = "steelblue4",
                          intgenes_size = 3,
                          labels_intgenes = TRUE,
                          labels_repel = F,
                          label_fontsize = 3){
  ma_df <- ma_df[ma_df$mean > 0, ]

  ma_df$logmean <- log10(ma_df$mean) # TO ALLOW FOR BRUSHING!!


  p <- ggplot(ma_df, aes_string(x = "logmean", y = "lfc", colour = "isDE"))

  if (!is.null(hlines)) {
    p <- p + geom_hline(aes(yintercept = hlines), col = "steelblue4", alpha = 1) +
      geom_hline(aes(yintercept = -hlines), col = "steelblue4", alpha = 1)
  }
  p <- p + geom_hline(aes(yintercept = 0), col = "red3", alpha = 1)

  p <- p + xlab(xlab) + ylab("log2 fold change")

  p <- p + geom_point(alpha = point_alpha)
  p <- p + scale_colour_manual(
    # name = paste0("FDR = ", FDR, "\n", "log2 fc = ", log2fc),
    name = "",
    values = c("TRUE" = sig_color, "FALSE" = "gray21"),
    labels = c("TRUE" = "DE", "FALSE" = "nonDE")
  )

  if (!is.null(ylim)) {
    p <- p + coord_cartesian(ylim = ylim)
  }

  if (!is.null(title)) {
    p <- p + ggtitle(title)+
      theme(plot.title = element_text(hjust = 0.5))
  }

  if (labels_intgenes && !is.null(intgenes)) {
      df_intgenes <- ma_df[ma_df$name %in% intgenes, ]
      df_intgenes$myids <- df_intgenes$name

    if(nrow(df_intgenes) > 0){
      p <- p + geom_point(data = df_intgenes, aes_string("logmean", "lfc"), color = intgenes_color, size = intgenes_size)

      if (labels_repel) {
        p <- p + ggrepel::geom_text_repel(
          data = df_intgenes, aes_string("logmean", "lfc", label = "myids"),
          color = intgenes_color, size = label_fontsize
        )
      } else {
        p <- p + geom_text(
          data = df_intgenes, aes_string("logmean", "lfc", label = "myids"),
          color = intgenes_color, size = label_fontsize, hjust = 0.25, vjust = -0.75
        )
      }
    }

  }

  if (add_rug) {
    p <- p + geom_rug(alpha = 0.3)
  }

  p <- p + theme_bw() + theme(plot.title = element_text(hjust = 0.5))
  p

}


#' Plot the fit normal for log2 fold change
#'
#' @param object A summarizedExperiment object
#' @param contrast The contrasts to plot fitting normal curve.
#'
#' @return
#' A histogram plot with a normal curve
#'
#' @export
#' @examples
#' # Load example
#' data(Silicosis_pg)
#' data <- Silicosis_pg
#' data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")
#'
#' # Differential test
#' ecols <- grep("LFQ.", colnames(data_unique))
#' se <- make_se_parse(data_unique, ecols,mode = "delim")
#' filt <- filter_se(se, thr = 0, fraction = 0.4, filter_formula = ~ Reverse != "+" & Potential.contaminant!="+")
#' norm <- normalize_vsn(filt)
#' imputed <- impute(norm, fun = "MinProb", q = 0.05)
#' diff <- test_diff(imputed, type = "control", control  = c("PBS"), fdr.type = "Storey's qvalue")
#' plot_norm_distribution(diff, contrast = "W4_vs_PBS")
#'
plot_norm_distribution <- function(object, contrast = get_contrast(object)) {
  row_data <- as.data.frame(rowData(object))
  cols_diff <- grep(paste("^",contrast, "_diff", sep = ""), colnames(row_data))


  x = row_data[, cols_diff]
  x = fun.outlier(x) %>% na.omit()
  fit <- fitnormal(x)
  Sigma = sqrt(fit$theta[2])
  mu = fit$theta[1]

  df <- as.data.frame(x)
  colnames(df) = gsub("_diff", "", colnames(row_data)[cols_diff])
  dat <- df %>% tibble::rownames_to_column() %>% tidyr::gather("contrast", "diff", -rowname)

  p <- ggplot(dat, aes(x = diff)) +
    geom_histogram(aes(y =..density..),
                   bins = 100,
                   colour ="black",
                   fill ="grey") +
    stat_function(fun = dnorm, args = list(mean = mu, sd = Sigma), color ="darkred") +
    # annotate("text", x = Inf, y = Inf, label = "Upper left", hjust = 2,
    # vjust = 2) +
    geom_text(data = data.frame(), aes(x = Inf, y = Inf,label = paste("Mean = ", round(mu, 2), "\nStd.Dev. =", round(Sigma, 2), sep = "")), hjust = 1.2, vjust = 1.2, size = 4) +
    labs(x = "log2 fold change") +
    theme_DEP2() +
    theme(legend.position = "none") +
    facet_wrap(~ contrast)
  return(p)
}


